
#include <list>
using namespace std;

#ifndef _RUNGEKUTTA4_H_
#define _RUNGEKUTTA4_H_

namespace Tuna {

template <typename T> class RungeKutta4;

template <typename T>
struct Typeinfo <RungeKutta4<T> > {
  typedef typename Typeinfo<T>::prec_t prec_t;
  enum { Dim = Typeinfo<T>::Dim };
};

//template<typename Interpolator_t, int Dim>
template<typename Interpolator_t>
class RungeKutta4
{
  typedef typename Typeinfo<Interpolator_t>::prec_t prec_t;
  typedef typename TunaArray<prec_t, Typeinfo<Interpolator_t>::Dim>::tiny floatTinyArray_t;
  typedef typename TunaArray<int, Typeinfo<Interpolator_t>::Dim>::tiny intTinyArray_t;
  typedef typename TunaArray<prec_t, Typeinfo<Interpolator_t>::Dim>::huge ScalarField_t;
///  Traits to define the type of arrays to be used.
///
///  floatTinyArray_t : Array of floats mostly used for deltas of meshes
///    
/*
    typedef typename TunaArray<double, Dim>::tiny floatTinyArray_t;
    typedef typename TunaArray<int, Dim>::tiny intTinyArray_t;    
    typedef typename TunaArray<double, Dim>::huge ScalarField_t;
*/
 public:
    RungeKutta4() { } ;
    ~RungeKutta4() { } ;
    
    inline floatTinyArray_t integrate(floatTinyArray_t &oldPos,
				      const floatTinyArray_t &deltas,
				      ScalarField_t &u,
				      ScalarField_t &v,
				      ScalarField_t &unext,
				      ScalarField_t &vnext,
				      prec_t dt,
				      floatTinyArray_t &vel_vort)
	{
	    UVW = i_object.interpolate(oldPos, deltas, u, v);
	    a(0) = (UVW(0) - vel_vort(0)) * dt;
	    a(1) = (UVW(1) - vel_vort(1)) * dt;

	    pos = oldPos;
	    pos(0) = pos(0) + a(0) * 0.5;
	    pos(1) = pos(1) + a(1) * 0.5;
	    UVW = i_object.interpolate(pos, deltas, u, v);
	    UVWnext = i_object.interpolate(pos, deltas, unext, vnext);
	    UVW(0) = ( UVW(0) + UVWnext(0) ) * 0.5;
	    UVW(1) = ( UVW(1) + UVWnext(1) ) * 0.5;
	    b(0) = (UVW(0) - vel_vort(0)) * dt;
	    b(1) = (UVW(1) - vel_vort(1)) * dt;
	    
	    pos = oldPos;
	    pos(0) = pos(0) + b(0) * 0.5;
	    pos(1) = pos(1) + b(1) * 0.5;
	    UVW = i_object.interpolate(pos, deltas, u, v);
	    UVWnext = i_object.interpolate(pos, deltas, unext, vnext);
	    UVW(0) = ( UVW(0) + UVWnext(0) ) * 0.5;
	    UVW(1) = ( UVW(1) + UVWnext(1) ) * 0.5;
	    c(0) = (UVW(0) - vel_vort(0)) * dt;
	    c(1) = (UVW(1) - vel_vort(1)) * dt;

	    pos = oldPos;
	    pos(0) = pos(0) + c(0);
	    pos(1) = pos(1) + c(1);
	    UVW = i_object.interpolate(pos, deltas, unext, vnext);
	    d(0) = (UVW(0) - vel_vort(0)) * dt;
	    d(1) = (UVW(1) - vel_vort(1)) * dt;
	    
	    newPos(0) = oldPos(0) + (a(0) + 2 * b(0) + 2 * c(0) + d(0)) / 6.0;
	    newPos(1) = oldPos(1) + (a(1) + 2 * b(1) + 2 * c(1) + d(1)) / 6.0;
	    
	    return newPos;
	}

    inline floatTinyArray_t integrate(floatTinyArray_t &oldPos,
				      const floatTinyArray_t &deltas,
				      ScalarField_t &u,
				      ScalarField_t &v,
				      ScalarField_t &w,
				      ScalarField_t &unext,
				      ScalarField_t &vnext,
				      ScalarField_t &wnext,
				      prec_t dt)
	{
	    UVW = i_object.interpolate(oldPos, deltas, u, v, w);
	    a(0) = UVW(0) * dt;
	    a(1) = UVW(1) * dt;
	    a(2) = UVW(2) * dt;

	    pos = oldPos;
	    pos(0) = pos(0) + a(0) * 0.5;
	    pos(1) = pos(1) + a(1) * 0.5;
	    pos(2) = pos(2) + a(2) * 0.5;
	    UVW = i_object.interpolate(pos, deltas, u, v, w);
	    UVWnext = i_object.interpolate(pos, deltas, unext, vnext, wnext);
	    UVW(0) = ( UVW(0) + UVWnext(0) ) * 0.5;
	    UVW(1) = ( UVW(1) + UVWnext(1) ) * 0.5;
	    UVW(2) = ( UVW(2) + UVWnext(2) ) * 0.5;
	    b(0) = UVW(0) * dt;
	    b(1) = UVW(1) * dt;
	    b(2) = UVW(2) * dt;
	    
	    pos = oldPos;
	    pos(0) = pos(0) + b(0) * 0.5;
	    pos(1) = pos(1) + b(1) * 0.5;
	    pos(2) = pos(2) + b(2) * 0.5;
	    UVW = i_object.interpolate(pos, deltas, u, v, w);
	    UVWnext = i_object.interpolate(pos, deltas, unext, vnext, wnext);
	    UVW(0) = ( UVW(0) + UVWnext(0) ) * 0.5;
	    UVW(1) = ( UVW(1) + UVWnext(1) ) * 0.5;
	    UVW(2) = ( UVW(2) + UVWnext(2) ) * 0.5;
	    c(0) = UVW(0) * dt;
	    c(1) = UVW(1) * dt;
	    c(2) = UVW(2) * dt;

	    pos = oldPos;
	    pos(0) = pos(0) + c(0);
	    pos(1) = pos(1) + c(1);
	    pos(2) = pos(2) + c(2);
	    UVW = i_object.interpolate(pos, deltas, unext, vnext, wnext);
	    d(0) = UVW(0) * dt;
	    d(1) = UVW(1) * dt;
	    d(2) = UVW(2) * dt;
	    
	    newPos(0) = oldPos(0) + (a(0) + 2 * b(0) + 2 * c(0) + d(0)) / 6.0;
	    newPos(1) = oldPos(1) + (a(1) + 2 * b(1) + 2 * c(1) + d(1)) / 6.0;
	    newPos(2) = oldPos(2) + (a(2) + 2 * b(2) + 2 * c(2) + d(2)) / 6.0;
	    
	    return newPos;
	}
    


    
 private:
    floatTinyArray_t UVW, UVWnext;
    floatTinyArray_t newPos, pos, a, b, c, d;
    Interpolator_t i_object;

};

} // Tuna namespace

#endif // _RUNGEKUTTA4_H_
